#include "MonteCarlo.h"

MonteCarlo::MonteCarlo(void){
}

MonteCarlo::~MonteCarlo(void){
}

double MonteCarlo::piNumber(int *idum, int begin, int end) {
	double x, y;
	int Pin = 0;
	int Ptotal = end - begin;
	int idum1, idum2;
	idum1 = *idum;
	long xix1, xix2,xix3;
	long yix1, yix2,yix3;
	double r1[98];
	double r2[98];
	int iff1, iff2;
	iff1 = iff2 = 0;
	for(int step=begin;step<=end;step++){
		x = rand1(&idum1, &xix1, &xix2, &xix3, r1, &iff1);
		y = rand1(&idum1, &xix1, &xix2, &xix3, r1, &iff1);
		if((pow(x,2)+pow(y,2))<=1)
			Pin++;
	}
	return (4.0*Pin/Ptotal);
}


double MonteCarlo::rand1(int *idum, long *ix1, long *ix2, long *ix3, double *r, int *iff){
	double temp;
	int j;
	if (*idum < 0 || iff == 0) {
		*iff = 1;
		*ix1 = (IC1-(*idum))%M1;
		*ix1 = (IA1*(*ix1)+IC1)%M1;
		*ix2 = (*ix1)%M2;
		*ix1 = (IA1*(*ix1)+IC1)%M1;
		*ix3 = (*ix1)%M3;
		for (j=1; j <=97; j++) {
			*(ix1) = (IA1*(*ix1)+IC1)%M1;
			*(ix2) = ((*ix1)*(*ix2)+IC2)%M2;
			r[j] = ((*ix1)+(*ix2)*RM2)*RM1;
		}
		*idum = 1;
	}
	*ix1 = (IA1*(*ix1)+IC1)%M1;
	*ix2 = (IA2*(*ix2)+IC2)%M2;
	*ix3 = (IA3*(*ix3)+IC3)%M3;
	j = 1 + ((97*(*ix3))/M3);
	if (j > 97 || j < 1)
		return 1;//erro

	temp=r[j];
	r[j] = ((*ix1)+(*ix2)*RM2)*RM1;
	return temp;
}
